Kaleidoscope of topological phases with multiple Majorana species 
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Exactly solvable lattice models for spins and non-interacting fermions provide fascinating examples of topolog- 
ical phases, some of them exhibiting the localized Majorana fermions that feature in proposals for topological 
quantum computing. The Chern invariant v is one important characterization of such phases. Here we look at 
the square-octagon variant of Kitaev's honeycomb model. It maps to spinful paired fermions and enjoys a rich 
phase diagram featuring distinct abelian and nonabelian phases with v — 0, ±1, ±2, ±3 and ±4. The v = ±1 
and v = ±3 phases all support localized Majorana modes and are examples of Ising and SU (2)2 anyon theories 
respectively. 



The study of topological quantum states has, in recent years, 
been boosted by proposals HH3) to use nonabelian anyons 
for topologically protected quantum computing. Nonabelian 
anyons in two-dimensional systems were first proposed in 
the context of strongly correlated systems, with the Moore- 
Read Fractional Quantum Hall state [4| being the most fa- 
mous example. Recent years have seen the advent of non- 
abelian anyons in effectively non-interacting systems. Pro- 
posals now exist for their realization in spinless chiral p-wave 
superconductors |5 1 and spin lattices |6 1, at interfaces between 
an s-wave superconductors and a topological insulator Q or 
half metal 0, and on 1-d semiconductor-superconductor het- 
erostructures I19 UTD . In all of the above examples the non- 
abelian quasi-particles are realized by isolated, topologically 
protected Majorana fermions. 

A necessary condition for isolated Majorana fermions is 
particle-hole symmetry of the type found in Bogoliubov-de 
Gennes (BdG) equations. A general understanding of where 
to find these conditions has been aided by the general sym- 
metry classification of non-interacting gapped bulk systems 
lfT2l4l5l . In particular, the symmetry class of 2d chiral p— and 
/—wave superconductors (class D in ifPfl ) is topologically 
non-trivial and can host topological bulk and edge-modes. 
The topological invariant characterizing the ground state in 
2d is the Chern number v, which maps to the set of all inte- 
gers Z. However,in the chiral p-wave case, for example, only 
Chern numbers v = (abelian phase) and v = ±1 (non- 
abelian phase) are uncovered from the full set Z. This is to 
be contrasted with the symmetry class of the Integer quantum 
Hall effect (class A), that lacks all symmetries. Here the Z 
invariant is fully realized (in theory) by filling some arbitrary 
number n of Landau level bands, which will be related to a 
Chern number v = ±n fl6l . 

One of our motivations for searching for higher Chern num- 
bers in systems with BdG symmetry is that phases with dif- 
ferent odd v host Majorana fermions of different topological 
properties. The Chern number dependence of the nonabelian 
braiding properties is discussed in Kitaev's seminal paper |6). 
A summary of the 16-fold way of bulk anyon models — so 
called because the braiding properties are determined by the 
value of v modulo 16 — is given in tables 1, 2 and 3 in that pa- 



per. The cases with odd v (but not those with even v) all sup- 
port localized Majorana modes at the vortices and thus have 
nonabelian exchange statistics [17]. We will see that the lat- 
tice models studied in this paper realizes 9 of these 16 bulk 
orders, including 4 of the cases with nonabelian anyons, in 
two parity-doublets v = ±1 and v — ±3. The only other 
example of Chern number v — ±3 that we are aware of oc- 
curs in the work of Mao and co-workers [18], which concerns 
an /-wave paired system of electrons moving in continuous 
spacetime. The Chern numbers v = ±1 and v = ±3 corre- 
spond to anyon models of the Ising and 5(7(2)2 type respec- 
tively. However, the edge behavior of the system is not fully 
determined by this classification and will depend on the ac- 
tual value of v (not its residue modulo 16). Nevertheless, by 
counting gapless edge modes for the various phases, we ar- 
gue that the v — ±1 edges are described by chiral Ising CFTs 
(Conformal Field Theories), while the v = ±3 edges support 
chiral SU{2) 2 CFTs. 

The model we study is the square-octagon (4-8) lattice 
model discussed by Yang et al. [19|. Like the triangle- 
dodecagon (3-12, "Yao-Kivelson") lattice model 11211 it is one 
of the possible generalizations in two dimensions JT9 1 of the 
original Kitaev honeycomb lattice model (6] l22T[30l . These 
exactly solvable 2d Ising-like spin models on trivalent lattices 
have turned out to be an ideal testing ground for topologi- 
cal properties predicted by the general symmetry arguments. 
The honeycomb model maps to spinless chiral p-wave system, 
and exhibits an abelian (v = 0) as well as nonabelian phases 
(v = ±1). Yang and coworkers observed [ 19 1 that the square- 
octagon variation has two abelian v = phases as well as the 
v = ±1 phases which open up near the critical boundary be- 
tween the abelian ones. In none of the mentioned models were 
ground state (vortex-free) sectors with higher Chern numbers 
encountered. However, since the square-octagon model has a 
larger unit cell than the model on the honeycomb lattice and 
can therefore be mapped to a fermionic paired system with 
an internal pseudo-spin degree of freedom, we suspected that 
it may support the rich variety of phases found in spin-triplet 
paired systems, see for example Ref.l20l 

We will see below that the the model supports an abun- 
dance of abelian and non-abelian phases with Chern num- 



bers 0, ±1, ±2, ±3, ±4 see Figure [3] Using the fermioniza- 
tion method introduced in Ref. 29 we have also checked di- 
rectly that this model can be mapped to a fermionic paired 
system with pairings of type s-, p-, d-, f- g- etc. along with 
spin-orbit couplings of corresponding winding numbers. In 
the absence of spin-orbit coupling — for example in spinless 
pairs — the correspondence between the Chern number and the 
angular momentum of the pairing is rather simple. A spinless 
chiral pairing of p-wave type corresponds in the weak-pairing 
phase to v = ±1, a chiral f-wave pairing instead to v = ±3 
etc. (In the strong-pairing phase one would in both cases find 
(i = 0.) The fact that we find phases with Chern numbers 
v = ±3 among others concurs generally with the f-wave pair- 
ing observed in Ref.[T8] However, we stress that in the model 
considered by us — which as mentioned contains intricate cou- 
plings of pseudospin and momentum with a variety of wind- 
ing numbers — we have not been able to discern any simple 
relation between dominant pairing and the Chern number. 

The outline of the paper is as follows. A general introduc- 
tion to the spin model will be accompanied by a description 
of the fermionic solutions, focusing in particular on the vortex 
free sector, which contains the ground state. We then describe 
the phase diagram and discuss how the edge spectrum, and in- 
deed the bulk spectrum at phase boundaries reflect the Chern 
number difference between topological phases. We also pro- 
vide details of how the contributions of each band in the model 
combine to give the total Chern number of each phase and de- 
scribe how these contributions change as the bands cross or 
touch. A general description of the edge spectra that occur 
between domains of the various topological phases (including 
vacuum) is then provided and we give an example where two 
domains with v = nevertheless have counter-propagating 
modes at the interface. We also show that these two abelian 
v = phases can be described with toric code like hamiltoni- 
ans, supported on different sub-lattice structures. 

The basic hamiltonian of the Kitaev-type trivalent spin 
model on the square-octagon (4-8) lattice is defined as 

h = -j z J2 aZ<jZ - J xJ2 ° X(jX ~ J vJ2 aV(jV (1) 

z-links x-links y-links 

where it should be understood that the z-links connect sep- 
arate squares of the lattice and the x and y-links within the 
squares are the horizontal and vertical slopes respectively, see 
Figure [T] We define a basic unit cell of the lattice around the 
four sites that make up a single square and label the sites in- 
side such a cell by n € 1, 2, 3, 4. Each spin site on the lattice 
can thus be specified using the position vector q of the unit 
cell and the index n. A unit cell and fundamental lattice vec- 
tors are shown in Figure [T] 

Within the model, as in honeycomb and 3-12 lattice mod- 
els, there are closed loop symmetries that we can gener- 
ate with overlapping link operators. The simplest of these 
are the octagonal and square plaquette opera- 

tors. These are defined pictorially in Figure [T] The fact that 
the hamiltonian commutes with all plaquette operators im- 
plies that we may choose energy eigenvectors | n) such that 
W^ a) = (n\W ( q a) \n) = ±1. For example, if = -1 
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FIG. 1: The square-octagon Kitaev lattice model. The horizontal 
links represent operators of the type a x a x , the vertical links repre- 
sent operators a y a y and the diagonal links represent the operators 
a z a z . Shown also are the fundamental unit cell with the sites inside 
it labeled by n 6 {1, ... , 4}, the lattice vectors n x and n y , and the 
plaquette operators W a . 

then we say that the state | n) carries a square vortex at q. By 
the vortex-sector we mean the subspace of the system Hilbert 
space with a particular configuration of vortices. The vortex- 
free sector for example is the subspace spanned by all eigen- 
vectors such that W^ 4) = 1 and = 1 for all q. The 

vortex-lattice sector, with Wq 4 ^ = — 1 and W q ^ = — 1 for all 
q, was shown in [3 1 1 to have a Fermi-surface and is thus a spin 
metal, like certain sectors of the 3-12 lattice variation |32|. 

In this type of system, it is possible to break time reversal 
symmetry explicitly, and still preserve the solvable nature of 
the system, by adding into the hamiltonian three site interac- 
tions of the form afa^a^. To find these terms we take over- 
lapping products of the adjacent links of the original hamilto- 
nian. These terms are shown in Figure [2] Crucial to obtain- 




FIG. 2: Time-reversal symmetry breaking terms. There are 8 octag- 
onal (O) terms and 4 square (Q) terms. We construct each term using 
the overlapping products of two adjacent o a a a terms. To fix the 
sign of each term, we use a clockwise ordering for the adjacent link 
operators. The resulting signs are noted inside the plaquettes. The O 
terms are added to the hamiltonian with an overall coupling constant 
2ki and the Q terms with an overall coupling K2 
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ing the rich phase structure shown in Figure [3] is the sign we 
place in front of each term. Following [19] we use the con- 
vention that the individual links are arranged in a clockwise 
order around each plaquette, which yields the signs shown in 
the figure. Two overall coupling constants k\ for the terms on 
square plaquettes and «2 for the terms on octagonal plaquettes 
allow us to vary the strength and nature of these T-breaking 
terms. 

To solve the model exactly one converts each vortex sec- 
tor of the model to a quadratic fermionic system. There are 
a number of different techniques that do this, each with their 
own distinct advantages. For the purposes of calculating the 
fermionic spectra it is easiest to use the Majorana represen- 
tation and conventions for this model used in 1191 1311 . It is 
important to note however that this methodology means we 
artificially enlarge the Hilbert-space and thus should project 
back to the physical subspace in order to calculate wavefunc- 
tions 0. Indeed, for the purpose of calculating wavefunctions 
it is arguably easier to use a Jordan- Wigner type fermioniza- 
tion procedure l27l l28l . We have elsewhere given details 
on a representation that uses a special Jordan-Wigner type 
string to resolve the fermionic vacuum as Toric code stabi- 
lized states [29 1 . Ground states calculated with this method- 
ology are BCS ground state wave functions over these Toric 
Code vacua and the full set of excited states can be described 
in terms of BdG quasi-particles over these ground states. 

We now focus on the vortex free sector, using the Majorana 
fermionic representation of ref. 19 After this fermionization, 
the Hamiltonian can be written as follows 



(2) 



where A\. = {a kA ,a k2 ,a k ^,a kA ) and the operators 
a k , n are the Fourier transformed Majorana fermion cre- 
ation/annihilation operators a kn = -j= ^ a q . n e lk ' q . The in- 
dices n label the four sites within a unit cell (see Figure [TJ, 
q is the center-of-mass coordinate of the unit cell, and a q ^ n is 
the operator which creates or annihilates a Majorana fermion 
at the site labeled by {q, n). The 4x4 matrix M can be de- 
composed as M = Mj + M K where 



Mj = 2i 
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The Chern invariant can be calculated as 



v=^.J d 2 kTr \P(k x ,k y ) 



dP dP dP dP 



(5) 



where P is the projection onto the negative energy eigenvec- 
tors of M. For multi-banded systems one can calculate a 
Chern integer for each band by using the same formula, but 
with P replaced by the projection to just that band. Each neg- 
ative energy band then contributes its own integer to the total 
Chern invariant. The Chern integers of the completely positive 
energy bands are just —1 times the analogues for the negative 
energy bands and therefore contain the same information. 

A change in the Chern number |5]l indicates a phase tran- 
sition (the reverse is not always true as Figure |4] illustrates). 
Changes in Chern number can occur when the positive and 
negative energy bands touch at one or several points. We can 
thus explore the phase diagram of the model by mapping out 
the Chern numbers and the gap as a function of the parameters 
{Jx, Jy, Jz, Ki, ^2). Clearly we can scale one of these param- 
eters to 1 (as long as it was non-zero to start with) and this 
means we have a 4-dimensional parameter space to deal with. 
We will present the full structure of the 4D phase diagram in 
a separate paper li33l . Here, we focus on a two-dimensional 
slice through parameter space, setting J z = 1, J = J x = J y 
and Ki — K2- In Fig.[3]we show the phases and Chern num- 
bers in this slice. 

The information in Figure [3] suggests that we can think of 
the topological phase transitions as an exchange of the Chern 
integers between positive and negative energy eigenmodes. 
Band touching within the negative energy bands also comes 
with a transfer of Chern integers, but here the total Chern in- 
variant is unchanged. In the model in question the two neg- 
ative energy bands touch whenever k = 0, k = ±J z /2 or 
k = ± J/2. These internal 'transitions' are also indicated in 
Figure [5] 

The dispersion relation for the energy E k:cky of the model 
can be obtained exactly, but is complicated and we will not 
write it down here. However, the critical lines shown in 
Figure [3] have simple analytical descriptions. We find that 
in all these cases the gap closes linearly and that the num- 
ber of isolated E kv .,k — points at a given transition dic- 
tates how the Chern invariant jumps between phases. For 
example, the critical lines of positive slope in the diagram, 
J = J (4k ± \ 2{J Z + 2k)) correspond to single zeros of the 
energy at {k x , k y ) = (0, tt). The same is true for the two lines 
of negative slope, J = |(— 4k ± \J2{J Z — 2k)) which have 
zero energy at the point {k x , k y ) = {tt, 0). 

The two hyperbolas in the diagram are given by k = 
±yJJl +2J 2 /2. In both cases, the bands touch when 
(k x ,k y ) — (0,0) and when (k Xl k y ) — (7r,7r) so that v 
jumps by 2 over these curves. The upper critical curve sur- 
rounding the v = ±3 phases and the lower critical curve 
surrounding the v — ±4 phases have the positive and neg- 
ative energy bands touch simultaneously at 4 points in mo- 
mentum space. These curves are given explicitly by k = 
±j(l + 2 J ± yf{l + 4J + 12J 2 )). The curve separating 
the v — — 4 phase from the v = phase and the curve 
separating the v = — 1 phase from the v = +3 phase have 
each four zero-modes of the form {tt, k y i), {tt, k y 2), {k x \, 0) 
and {k X 2,Q) with k x \,k x i,k v \,k y i ^ {0,7r}. Similarly, the 
curves separating the v = 4 phase from the v = phase and 
the v = +1 phase from the v = — 3 phase have -Efc x ,fc = at 
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FIG. 3: Phase diagram of the square-octagon model. In this figure J z = 1, J = J x = J v and k = Ki = ft2- The dark lines indicate phase 
transitions where the higher negative energy band and the lower positive energy band touch. Inside each phase the Chern numbers of the two 
negative bands are also indicated. The top number indicates the higher negative energy band, the one lying closest to E — 0. The dotted lines 
indicate where these negative energy bands touch. In this case the internal transfer of Chern integers does not change the total Chern invariant. 
All phases with total Chern number ±3, ±4 have topological orders for which we do not know of previous lattice realizations. 



(0,k yl ),(0,k y2 ),(k xl ,ir) and {k x2 ,n). 

As proved in [6|, each vortex in the odd Chern number 
phases is accompanied by a single localized zero energy Ma- 
jorana fermion excitation (assuming the vortices are well- 
separated). These localized Majorana fermions in the dif- 
ferent phases with odd v make the vortices into nonabelian 
anyons and we obtain here 4 out of the 8 different types non- 
abelian anyons which occur in Kitaev's 16-fold way for free 
fermionic bulk topological orders (in 2 pairs of types related 
by parity). Specifically, the v = ±1 modes correspond to the 
so called Ising anyon theory and the v — ±3 modes corre- 
spond to the SU (2)2 anyon theory. These zero energy excita- 
tions give rise to a ground state degeneracy and allow for the 
possibility of nonabelian statistics and fault tolerant quantum 
computation QJj3). As far as we know, this is the first time 
that different types of nonabelian anyons carrying Majorana 
modes have been found in a lattice model. We find zero en- 
ergy modes at vortex excitations also in some of the abelian 
phases of the model, whose topological orders cover 5 of the 
8 abelian possibilities. In these cases however, the modes are 
full Dirac modes (consisting of two Majorana modes) and are 
thus not topologically protected from local error processes. 
We will discuss one of these Dirac modes briefly below. 

The multitude of topological phases in the model allow one 
to set up a variety of interesting domain wall scenarios. The 
general result of Hatsugai [34| always applies here: The net 
number of left movers minus the number of right movers is 
equal to the difference in Chern number of the phases on ei- 
ther side of the boundary. In the simplest cases, for example 
for open boundaries, where the topological medium just ends, 
we find that this is usually satisfied in the simplest possible 
way; the edge theory is chiral and has a number of zero modes 
equal to \v\. This supports the idea that the open edge theo- 



ries of the v — ±1 and v = ±3 phases are Ising and SU(2)2 
CFTs, the simplest CFTs which satisfy the bulk-boundary cor- 
respondence. 

For edge theories on domain walls the same is often true but 
not always. More generally we find that the number and type 
of zero modes that exist on a domain wall always reflects the 
critical bulk theory that exists on the phase boundary between 
the phases on either side of the wall. If one examines the edge 
spectra in a system with one phase domain below another in 
such a way that the translational invariance in the x direction 
is maintained, one finds that the wall contains a band of states 
that intersect the E = line at value of k x that are consistent 
with the k x values of the zero modes that occur at the bulk 
phase transition. If the two phases are separated in the phase 
diagram by more than one phase boundary, then the edge will 
reflect all the phase boundaries crossed by an interpolating 
path. Examples of this behavior are given in Figures [4] and [5] 

Careful inspection of the phase diagram reveals a number 
of situations where phases with the same Chern number touch 
at the intersections and touchings of critical lines. In these 
cases, we can have domain wall edge states between phases 
at the same Chern number. Such domain wall edge states are 
simply inherited from the boundaries between phases with dif- 
ferent Chern numbers that meet at these multicritical points 
and hence they are stable against perturbations which do not 
remove the multicritical points. However, some of the mul- 
ticritical points and corresponding edge states can be easily 
removed by local perturbations to the Hamiltonian, while oth- 
ers are topologically protected, that is, stable against any local 
perturbation which preserves particle-hole symmetry. We now 
give an example of each case. Following Yang and cowork- 
ers |[19| we call the phase around the n = 0,\J\ < J z /\/2 
domain the A z phase and around the n = 0, | J| > J 2 /\/2 
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FIG. 4: The central figure shows an excerpt of the central region 
of the phase diagram in Fig. [5] The domain wall edge spectra be- 
tween two phases reflect all the bulk phase transitions that occur in 
the phase diagram along a continuous line connecting the two phases. 
Low energy edge modes associated with the red edge shown in red 
(in green for the green edge). We show here the bulk phase diagram 
at a A xy — B and a A z — B transition and how it is reflected in the 
domain wall edge spectra. One can also see here how the representa- 
tive edge spectra from A xy — B edge and A z — B edges are in some 
sense combined at an A xy — A z edge. The A xy and A z phases can 
both be perturbatively mapped to Toric Code type systems that are 
supported on different sublattices. 



line the A xy phase. It turns out that we can make the tetra- 
critical points that occur where the critical lines of positive 
slope cross and where the critical lines of negative slope cross 
vanish by varying J y away from J x , in the process connect- 
ing the v — phases on the left and right to the A xy phase. 
In contrast, the tetra-critical points between the A xy and A z 
phases are robust features that cannot be so easily removed. 
In Figure[4]we show the edge spectrum for two horizontal do- 
main walls between A z phase and the upper A xy phase. We 
see that on a single edge we have a single band of both left 
and right moving states that intersects E — at k x =0 and 
k x — ir. We now argue for the stability of these edge bands. 

First of all, if the system is wide enough in the y-direction, 
the edge modes on opposite edges in Fig. |4] decouple expo- 
nentially and tunneling between the edges cannot lift the zero 
modes. Backscattering between the leftmoving and rightmov- 
ing zero modes on each edge could in principle lift these 
modes, but this would violate particle-hole symmetry. This 
is easiest to see for perturbations which preserve translational 
symmetry. In this case, because we have E(k) = —E(—k) 



the modes at k x — and k x = ir must remain at zero energy. 
If translational symmetry is broken (for example by disorder), 
a slightly more delicate argument is needed. In this case k x is 
no longer a good quantum number, but if we imagine that the 
disorder was introduced adiabatically, we may still label the 
energy eigenstates by the former fc^-eigenvalues. The disor- 
der may now make the energy of the k x — and k x = ir states 
nonzero, but in order to preserve particle hole symmetry, it 
must still preserve the property that every state at positive en- 
ergy has a partner at negative energy. Therefore, in the con- 
tinuum limit, where k x takes continuous values and energy is 
a continuous function of k x , there will always be zero modes 
on the domain wall, though not necessarily at k x — and 
k x = ir. In fact, there must always be a nonzero even number 
of zero modes, since zero modes can adiabatically appear and 
disappear in even numbers only. 

Of course k. x is discrete for any finite system size and these 
arguments do not guarantee the existence of modes with ex- 
actly zero energy, just increasingly soft modes for growing 
system size. This is in fact how the system behaves when there 
is a finite and even number of unit cells in the x-direction. 
However, if there is an odd number of unit cells in the re- 
direction, the zero energy modes at k x = and k x = ir 
actually belong to different topological sectors and therefore 
backscattering between the k x = and k x = it modes would 
violate conservation of a topologically conserved quantum 
number. In this situation we can have a robust Majorana zero 
mode forming on an interface between v = phases even at 
finite system size. 

It is gratifying that both Abelian phases in the model can be 
analysed using degenerate perturbation theory. For the v = 
phase found in the center of the figure (the isolated z-link limit 
( J z ~> (J x + Jy ) 1 ^ 2 ) it was shown in [19] that the low-energy 
effective system, at the 4th order of perturbation theory, is that 
of an abelian toric code on a square lattice, 



c 72 t2 

0J x J y \ " r»r(8) 



To the 4th order the constant terms can be calculated to be 
E = N{\J Z \ + ^ + 3^3) where N is number of z- 
links or fermions and is thus twice the number of unit cells. 
In addition, we have found that such a mapping is possible 
also in the isolated square limit (J z < (Jj + J2)V2). Each 
isolated square is a miniature spin chain with alternating X 
and Y links as given by eq. [T] In the absence of a vortex 
on the plaquette, it has a doubly degenerate ground state. This 
degenerate subspace acts as an effective spin in much the same 
way as the ground states of the isolated z-dimers above do. 
Perturbing away from the fully isolated regime, we can derive, 
at the 4th order, the effective hamiltonian 



(8) 



(6) 



where r = 2{J% + J^) 1 / 2 and E 



(8) 

suitable basis the operators W q are of the form afa^a^a^. 



iV(r + i + ^). In a 
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Here the subscript labels indicate left, right, up and down re- 
spectively and the operators a operate on the low energy anti- 
periodic sector of the isolated 4-spin chain, see Fig. [4] In this 
regime, excitations on the octagonal plaquettes can be inter- 
preted as electric (e) and magnetic (to) particles, in a checker- 
board pattern as in ref. (6). We see then that we have a toric 
code type phase, but on a different lattice than that found for 
the isolated z-link limit 

Further distinctions between these Abelian phases can be 
made if we examine the system's square vortex excitations. 
Note first that the spectrum of the 4-spin with a vortex exci- 
tation is four fold degenerate. This degeneracy is not lifted 
when we perturb away from the isolated square limit and so 
we can infer that in this A xy phase there exists a single mass- 
less Dirac fermion mode in the presence of a square-vortex 
excitation. It is important to note that this degeneracy is de- 
pendent on the fact that the chain is in an exact eigenstate of 
the square plaquette operator and hence it is not topologically 
protected (it can be broken by local perturbations). 
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FIG. 5: The edge spectra (a) between the C (u = +3) phase and the 
A z (v = 0) phase, (b) between the C (y = +3) phase and the A xy 
(i> = 0) phase and (c) between the C (y = +3) phase and the B 
(y = — 1) phase 

The domain wall between the abelian A z and A xy phases 
at v = is not the only one which shows counter-propagating 
zero modes. To illustrate this fact, we include figure [5] which 
shows that domain walls involving the nonabelian C phase, 
at v = 3, exhibit similar phenomena. The C — A z (G — E) 
boundary spectra shown behave as would be naively expected, 
with 3 (4) zero modes of the same chirality on each of the 
respective edges (shown in red and green). However, the C — 
A xy boundary actually retains the zero modes at k x = and 
k x = 7r that comes from either taking the path C — B — A xy 



or C — A z — A xy through the phase diagram. 

In conclusion, in this paper we have seen that the square- 
octagon lattice variation of the Kitaev honeycomb model has 
a startling richness of topological phases. In particular it has 
two essentially different types of nonabelian anyons at Chern 
numbers v = ±1 and v = ±3 and it realizes 9 of the 16 
bulk topological orders allowed by Kitaev's 16-fold way. We 
demonstrated how the bulk phase diagram of the system can 
be used to predict the behavior of the edge spectra between 
different topological domains. We encounter various situa- 
tions where phases with the same Chern number have a non- 
zero (but even) number of low energy edge bands due to the 
interface. An examination of the perturbation theory of the 
abelian phases reveals that distinct abelian phases can both 
be mapped to Toric Code type systems that are supported on 
different sub-lattice structures. 

The perhaps most interesting aspect of the model is that it 
gives an example of a system with Majorana modes that seizes 
the possibility of higher Chern numbers allowed from general 
symmetry classifications of Altland and Zirnbauer ifTZl [T3l . 
Compared to the Kitaev honeycomb model that maps to spin- 
less paired fermions, the pseudo fermion spin in the square- 
octagon model doubles the number of bands and from this 
simple doubling of bands we could expect Chern numbers 
from up to ±2. The real surprise is the occurrence of 
v = ±3, ±4 which are beyond the expected doubling. In 
particular, we have found instances where a single band can 
contribute with v = ±3. 

A deeper understanding of how in detail the multi-band sit- 
uation allows for increased Chern numbers is left for future 
study, but it may be imagined that, given a sufficient num- 
ber of sites in the unit cell of the lattice, it should be possible 
to realize all 16 possible bulk topological orders for gapped 
Majorana systems within the same trivalent spin models. In 
this respect, some progress has also been made using Bloch 
Hamiltonians derived from the vortex sectors of the Kitaev 
honeycomb model and related continuum systems, see for ex- 
ample Ref. [35 1 . It is also interesting to speculate on the sit- 
uation for gapped phases in 3d variations [36] of Kitaev type 
models. The winding number is in this case given by a Wess- 
Zumino-Witten term[ 14 1. As the number of bands is typically 
higher, at least four, it seems probable that winding numbers 
other than v = 0,±1 could be found. The general symme- 
try arguments lfT4l provide a situation that is almost inverse 
to that in 2d. Only one of the classes with a Z bulk comes 
with broken T-symmetry, and this class — AIII — is probably 
less relevant to the Kitaev type models as it lacks particle- 
hole symmetry. The other two classes with Z bulks, namely 
Dili and CI, require particle-hole symmetry and time-reversal 
symmetry. Nevertheless, it remains to be seen whether higher 
winding numbers can actually be realized, possibly by extend- 
ing these T-invariant models with, for example, additional 4- 
spin interactions. 
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